#include "PowerSpectrum.h"

#include <iostream>
#include <fstream>
using std::endl;
using std::cout;

#include "dbg.h"

std::ostream* dbgout=0;
bool XDEBUG=false;

int main(int argc, char* argv[])
{
  dbgout = &std::cout;

  // Om_m  Om_DE  Om_b  h  A_s  n  w  wa
  const Universe u(0.3,0.7,0.024,0.7,0.8,1.0,-1.0,0.0);

  double meanz = 0., meanw = 0., sumpbz = 0., sumpb = 0.;
  cout<<"#    z    w    p(z)   p(w)\n";
  for(double z = 0.05; z < 3.; z+=0.05) {
    double w = u.w_z(z);
    double pb = u.pb(w);
    double pbz = pb*u.dwdz(z);
    cout<<z<<"  "<<w<<"  "<<pbz<<"  "<<pb<<endl;
    meanz += pbz*z;
    meanw += pb*w;
    sumpbz += pbz;
    sumpb += pb;
  }

  meanz /= sumpbz;
  meanw /= sumpb;

  cout<<"Mean z = "<<meanz<<endl;
  cout<<"Mean w = "<<meanw<<endl;

  return 0;
}
